rm(list=ls())
gc()
require(foreign)
require(ineq)
require(dplyr)
library(Hmisc)
require(corrplot)
require(systemfit)
require(cem)
require(BayesTree)
require(stargazer)
require(texreg)
require(lavaan)

setwd(paste(substr(getwd(),1,regexpr("Dropbox",getwd())[1]-1),"Dropbox/coauthors/AMT_Bisbee_Larson/bisbee_larson_replication_final",sep=""))
getwd()
source("./tools/leg2.R")
##############################################################################################
##
##  File Name:    rels_results.R
##
##  Input Files:  rels_donate_cem.dta
##                final_data_long.dta
##  Output Files: rels_donate_cem.pdf             (Figure 8)
##                SI_sem_donate.pdf               (SI Figure 4)
##
##  Purpose: This file summarizes the relationship results for donation behavior in 
##            rels_donate_cem.pdf (Figure 8 in the text) as well as applying structural
##            equation modeling (SEM) to the data to generate SI_sem_donate.pdf (Figure 4
##            in the supporting information).
##
##############################################################################################

betas <- read.dta("./relationships/rels_donate_cem.dta")

###############################################################################################
# rels_donate_cem.pdf: Figure 8
###############################################################################################
betas <- betas[complete.cases(betas),]
betas$diffs <- betas$betaon - betas$betaoff
betas$ses <- sqrt(betas$seon^2+betas$seoff^2)
betas$tstats <- betas$diffs/betas$ses
coefs <- as.data.frame(cbind(as.numeric(betas$betaon),as.numeric(betas$betaoff)))
coefs$diff <- abs(as.numeric(coefs[,1])-as.numeric(coefs[,2]))
ses <- as.data.frame(cbind(betas$seon,betas$seoff))
y.axis <- length(betas$betaon):1

order1 <- order(betas$diffs)
coefs1 <- coefs[order1,]
ses1 <- ses[order1,]
labels1 <- betas[order1,"names"]
betas$tie <- sapply(1:length(betas[,"names"]),function(i) strsplit(betas[,"names"]," ")[[i]][length(strsplit(betas[,"names"]," ")[[i]])])
betas$tie[which(betas$tie=="Weakes")] <- "Weakest"

betas$names2 <- sapply(1:nrow(betas),function(i) substr(betas$names[i],1,stop = gregexpr(" ",betas$names[i])[[1]][length(gregexpr(" ",betas$names[i])[[1]])]-1))

betas1 <- betas %>% select(names2,tie,diffs,ses)
betas1 <- reshape(betas1,timevar="tie",idvar=c("names2"),direction="wide")
xlim1 <- c(min(betas$diffs-2*betas$ses),max(betas$diffs+2*betas$ses))

onmin <- coefs[,1]-2*ses[,1]
offmin <- coefs[,2]-2*ses[,2]
minmin <- NA
for(i in 1:length(onmin)) {
  minmin[i] <- min(onmin[i],offmin[i])
}

order2 <- order(minmin)
coefs2 <- coefs[order2,]
ses2 <- ses[order2,]
labels2 <- betas[order2,"names"]
ord <- nrow(betas1):1

# Sprucing up the figure
betas2 <- betas %>% select(betaon,seon,betaoff,seoff,names2,tie,tstats)
betas2 <- reshape(betas2,timevar="tie",idvar=c("names2"),direction="wide")
#betas2[which(betas2$names2 == "Preferred Interaction" | betas2$names2 == "Tie Strength"),which(grepl("beta",colnames(betas2)))] <- betas2 %>% filter(names2 == "Preferred Interaction" | names2 == "Tie Strength") %>% select(which(grepl("beta",colnames(betas2))))*-1
betas2[,which(grepl("tstat",colnames(betas2)))] <- abs(betas2[,which(grepl("tstat",colnames(betas2)))])
y.axis <- 1:length(betas2[,1])
betas2$names2 <- c("Job Search","Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)","Education Hom.","Religious Hom.","Political Hom.","Class Hom.","Tie Strength","Pers. Crisis","Pers. Success","Prof. Crisis","Prof. Success","Pref. Interaction")
order <- match(c("Political Hom.","Religious Hom.","Class Hom.","Education Hom.",
           "Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)",
           "Pers. Success","Pers. Crisis","Pref. Interaction","Tie Strength",
           "Job Search","Prof. Success","Prof. Crisis"),betas2$names2)

betas2 <- betas2[rev(order),]
pdf("./1_Figures/rels_donate_cem.pdf",height=5,width=6)
nf <- layout(matrix(c(1,2,3,4,5,6,7,7,7,7,7,7),2,6,byrow=T),widths=c(2,1.4,1.4,1.4,1.4,1.4),heights=c(7,.75))
#layout.show(nf)
par(mar=c(2,5,3,0))
plot(rep(0,length(y.axis)), y.axis, type = "n", axes = F, xlab = "", ylab = "",ylim = c(.5,length(y.axis)+.5), main = "")
axis(2, at = y.axis, label = betas2[,1], las = 1, tick = FALSE, cex.axis =1.3,pos = .9)

par(mar=c(2,.5,3,.5))
for(i in c("Strongest","Strong","Weak","Weakest","Clique")) {
  min <- min(betas2[,paste("betaon.",i,sep="")]-3*betas2[,paste("seon.",i,sep="")],
             betas2[,paste("betaoff.",i,sep="")]-3*betas2[,paste("seoff.",i,sep="")])
  max <- max(betas2[,paste("betaon.",i,sep="")]+7*betas2[,paste("seon.",i,sep="")],
             betas2[,paste("betaoff.",i,sep="")]+7*betas2[,paste("seoff.",i,sep="")])
  plot(betas2[,paste("betaon.",i,sep="")],y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.2, 
       xlim=c(min,max),ylim = c(.5,length(y.axis)+.5), main = i)
  axis(1, at = seq(min,max,(max-min)/10), 
       labels = sapply(0:10, function(x) round(min+x*((max-min)/10),1)),tick = T,cex.axis = .75, mgp = c(2,.7,0))
  segments(betas2[,paste("betaon.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis, 
           betas2[,paste("betaon.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis, lwd =  1,col=rgb(0,0,0,alpha=150,maxColorValue=255))
  segments(betas2[,paste("betaon.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis-.05, 
           betas2[,paste("betaon.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=150,maxColorValue=255))
  segments(betas2[,paste("betaon.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis-.05, 
           betas2[,paste("betaon.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seon.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=150,maxColorValue=255))
  
  segments(betas2[,paste("betaoff.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis, 
           betas2[,paste("betaoff.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis, lwd =  1,col=rgb(0,0,0,alpha=50,maxColorValue=255))
  segments(betas2[,paste("betaoff.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis-.05, 
           betas2[,paste("betaoff.",i,sep="")]-qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=50,maxColorValue=255))
  segments(betas2[,paste("betaoff.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis-.05, 
           betas2[,paste("betaoff.",i,sep="")]+qnorm(1-0.025/length(y.axis))*betas2[,paste("seoff.",i,sep="")], y.axis+.05, lwd =  1,col=rgb(0,0,0,alpha=50,maxColorValue=255))
  points(betas2[,paste("betaon.",i,sep="")],y.axis,pch=21,bg=rgb(0,0,0,alpha=150,maxColorValue=255),col=rgb(0,0,0,alpha=150,maxColorValue=255))
  points(betas2[,paste("betaoff.",i,sep="")],y.axis,pch=21,bg=rgb(0,0,0,alpha=50,maxColorValue=255),col=rgb(0,0,0,alpha=50,maxColorValue=255))
  stars <- ifelse(round(2*pt(-abs(betas2[,paste("tstats.",i,sep="")]),df = 450),2) < .05,
                  paste(round(2*pt(-abs(betas2[,paste("tstats.",i,sep="")]),df = 450),2),"*",sep=""),
                  round(2*pt(-abs(betas2[,paste("tstats.",i,sep="")]),df = 450),2))
  text(max,y.axis,stars,cex=.8,adj = .9)
  abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing
}

#add legend (manually) to identify which dots denote model 1 and which denote model 2
par(mar=c(0,0,0,0))
plot(rep(0,length(y.axis)), y.axis, type = "n", axes = F, xlab = "", ylab = "",ylim = c(.5,length(y.axis)+.5), main = "")
legend("center", c("Offline","Online","*  Sig @ 95%","**  Sig @ MC"), 
       lwd=c(1,1,NA,NA),horiz=T,pch=c(21,21,NA,NA),col=c("grey60","black",NA,NA),
       pt.bg = c(rgb(0,0,0,alpha=50,maxColorValue=255),rgb(0,0,0,alpha=250,maxColorValue=255),"black"),bg="white",bty='o', cex=.9)
dev.off()




#################################################################################################
##  SI_sem_donate.pdf: Supporting Information Figure 4
#################################################################################################
dat <- read.dta("./tools/final_data_clean.dta")

dims <- c("jobsearch","contribute","garner","personalgain","grpgain","homo_pol","homo_rel","homo_educ",
          "homo_ls","pers_cris_synth","pers_succ_synth","prof_cris_synth","prof_succ_synth",
          "seat_synth2_lst","seat_synth2_mst","interaction")

coefs <- ses <- matrix(NA,nrow = length(dims),ncol = 5)
i = 1
for(tie in c("clique","weakest","weak","strong","strongest")) {
  model <- paste("std_donate_",tie," ~ ",paste(paste("std_",dims,"_",tie,sep=""),collapse=" + "),sep = "")
  fit <- sem(model,data = dat)
  coefs[,i] <- fit@ParTable$est[which(grepl("^~$",fit@ParTable$op))]
  ses[,i] <- fit@ParTable$se[which(grepl("^~$",fit@ParTable$op))]
  i = i+1
}

rownames(coefs) <- rownames(ses) <- dims
colnames(coefs) <- colnames(ses) <- c("clique","weakest","weak","strong","strongest")

dat.long <- read.dta("./tools/final_data_long.dta")
dims <- c("jobsearch","contribute","garner","personalgain","grpgain","homo_pol","homo_rel","homo_educ",
          "homo_ls","pers_cris_synth","pers_succ_synth","prof_cris_synth","prof_succ_synth",
          "seat_synth2_lst","seat_synth2_mst","interaction")
model.long <- paste("std_donate_ ~ ",paste(paste("std_",dims,"_",sep=""),collapse=" + "),sep = "")
fit.long <- sem(model.long,data = dat.long)
coefs.long <- fit.long@ParTable$est[which(grepl("^~$",fit.long@ParTable$op))]
ses.long <- fit.long@ParTable$se[which(grepl("^~$",fit.long@ParTable$op))]

ords <- order(coefs.long)
labs <- c("Job Search","Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)","Political Hom.",
          "Religious Hom.","Education Hom.","Class Hom.","Pers. Crisis","Pers. Success",
          "Prof. Crisis","Prof. Success","Least in Common","Most in Common","Pref. Interaction")

pdf("./1_Figures/SI_sem_donate.pdf")
nf <- layout(matrix(c(1,2,3,4,5,6,7),1,7,byrow=T),widths=c(2,3,1))
#layout.show(nf)
par(mar=c(2,5,3,0))
plot(rep(0,length(coefs.long)), 1:length(coefs.long), 
     type = "n", axes = F, xlab = "", ylab = "",main = "")
axis(2, at = 1:length(coefs.long), label = labs[ords], las = 1, tick = FALSE, cex.axis =.9,pos = .9)

par(mar=c(2,0,3,1))
pchs <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,19,21)
cexs <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,.6,1.2)
cols <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,rgb(0,0,0,.1),rgb(0,0,0,.8))
ltys <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,1,1)
bgs <- ifelse(abs(coefs.long[ords]/ses.long[ords]) < 2.33,NA,"white")
plot(coefs.long,1:length(coefs.long),type='n',xlab="Coefficient",ylab="",yaxt='n',xlim=c(-.1,.25),main="Full Data")
segments(coefs.long[ords] - 2*ses.long[ords],1:length(coefs.long),coefs.long[ords] + 2*ses.long[ords],col=cols,lty=ltys)
points(coefs.long[ords],1:length(coefs.long),pch=pchs,col=cols,cex=cexs,bg=bgs)
abline(v=0,lty=2)
par(mar=c(2,0,3,.05))
ties <- c("Clique","Weakest","Weak","Strong","Strongest")
# plot(rep(0,length(coefs.long)), 1:length(coefs.long),xlim=c(-.2,.4),
#      type = "n", xlab = "", yaxt='n',ylab = "",main = "By Tie")
yadj <- 0
for(x in 1:5) {
  pchs <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,19,21)
  cexs <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,.6,1.2)
  cols <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,rgb(0,0,0,.1),rgb(0,0,0,.8))
  ltys <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,1,1)
  bgs <- ifelse(abs(coefs[ords,x]/ses[ords,x]) < 2,NA,"white")
  plot(rep(0,length(coefs.long)), 1:length(coefs.long),xlim=c(-.25,.35),
       type = "n", xlab = "", yaxt='n',ylab = "",main = ties[x],xaxt='n',bty='n')
  segments(coefs[ords,x] - 2*ses[ords,x],(1:length(coefs[,1]))+yadj,coefs[ords,x] + 2*ses[ords,x],col=cols,lty=ltys)
  points(coefs[ords,x],(1:length(coefs[,1]))+yadj,pch=pchs,cex=cexs,col = cols,bg=bgs)
  #yadj = yadj - .15
  abline(v=0,lty=2)
  abline(v = .37,col=ifelse(x == 5,"black",rgb(0,0,0,.3)))
  abline(v = -.27,col=ifelse(x == 1,"black",rgb(0,0,0,.3)))
  abline(h = 16.6)
  abline(h = .4)
}
dev.off()
